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Three-dimensional  microstructure  of  mixed  ionic  and  electronic  conducting  cathode, 
Lao.6Sr0.4Coo.2Feo.803_5  (LSCF6428),  is  obtained  by  a  dual-beam  focused  ion  beam-scanning  elec¬ 
tron  microscope,  and  its  overpotential  is  predicted  by  the  lattice  Boltzmann  method.  Gaseous,  ionic 
and  electronic  transport  equations  coupled  with  electrochemical  reaction  at  the  gas/solid  interface  in 
the  three-dimensional  microstructure  are  solved  with  an  assumption  of  local  equilibrium  in  the  solid 
oxide.  The  gas  transport  is  modeled  by  the  dusty  gas  model.  The  numerical  simulation  is  conducted 
under  the  current  density  conditions  of  0.01,  0.05,  0.1  and  0.2  A/cm2.  Predicted  cathode  overpotentials 
agreed  well  with  the  experimental  results.  However,  predicted  overpotential  was  very  large  at  O2  =  20%, 
T=973  K  and  i  =  0.2  A/cm2  case  due  to  the  decline  of  ionic  conductivity  at  low  oxygen  partial  pressure. 
Three-dimensional  chemical  potential  and  current  vector  distributions  inside  LSCF  microstructure  are 
presented.  Ionic  and  electronic  current  stream  lines  are  uniform  and  smooth,  which  indicates  good  ionic 
and  electronic  conductions  as  well  as  wide  electrochemically  active  areas  inside  the  LSCF  microstructure. 
Present  method  will  be  an  effective  tool  for  investigating  local  oxygen  potential  field  which  affects  local 
reactions,  diffusions  and  physical  properties  of  the  MIEC  cathodes. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

A  solid  oxide  fuel  cell  (SOFC)  is  expected  to  be  a  promising 
energy  conversion  device  in  the  future  because  of  its  high  efficiency 
and  fuel  flexibility  [1  ].  Cost  reduction  and  long-term  durability  are 
the  major  issues  for  the  development  of  SOFC  system.  One  of  the 
promising  approaches  is  the  reduction  of  operating  temperature, 
which  promotes  utilization  of  low  cost  metals  in  the  system.  On  the 
other  hand,  low  temperature  operation  leads  to  poor  performance 
of  the  cell  due  to  the  degradations  of  conductivities  and  polarization 
characteristics.  Thus,  great  efforts  have  been  paid  for  develop¬ 
ing  efficient  and  stable  materials  for  the  electrodes.  Mixed  ionic 
and  electronic  conductors  (MIECs)  such  as  Lai_xSrxCo03_5  (LSC) 
or  La1_xSrxCoi_yFey03_5  (LSCF)  are  expected  to  be  very  effective 
in  reducing  polarization  losses  of  the  electrodes  [2,3].  Polarization 
characteristics  of  MIEC  electrodes  depend  on  electrochemical  reac¬ 
tion  rates  at  the  gas/solid  interface,  and  diffusions  of  ion,  electron 
and  gas  species,  which  are  strongly  affected  by  the  microstructure. 
Recently,  direct  measurements  of  three-dimensional  SOFC  elec¬ 
trode  microstructures  have  been  conducted  by  focused  ion  beam 
scanning  electron  microscopy  (FIB-SEM)  [4-7]  and  X-ray  computed 
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tomography  (XCT)  [8,9].  Through  these  three-dimensional  mea¬ 
surements,  important  microstructural  parameters  such  as  three 
phase  boundary  (TPB)  length  and  tortuosity  factors  can  be  obtained. 
For  example,  in  Gostovic  et  al.  [5],  three-dimensional  microstruc¬ 
ture  of  LSCF  cathode  was  reconstructed  by  FIB-SEM  and  the 
relationship  between  cathodic  polarization  and  surface  area  or  tor¬ 
tuosity  has  been  investigated. 

In  order  to  discuss  quantitative  relationship  between  electrode 
microstructures  and  their  polarization  characteristics,  numerical 
simulations  are  expected  to  provide  useful  information  which  can¬ 
not  be  obtained  from  experiments.  Joshi  et  al.  [10]  performed  a 
multi-component  lattice  Boltzmann  method  (LBM)  simulation  in 
a  two-dimensional  porous  media.  Asinari  et  al.  [1 1  ]  also  used  LBM 
for  solving  H2  and  H20  diffusion  inside  the  micro  pores.  Suzue  et  al. 

[12]  conducted  a  three-dimensional  LBM  simulation  which  solves 
the  species  transport  coupled  with  the  electrochemical  reaction  in  a 
stochastically  reconstructed  anode  microstructure.  Shikazono  et  al. 

[13]  used  Ni/YSZ  anode  microstructure  reconstructed  by  FIB-SEM 
and  solved  transport  equations  and  electrochemical  reactions  at 
TPB  using  LBM.  Three-dimensional  potential  and  current  vector  dis¬ 
tributions  inside  actual  Ni/YSZ  anode  are  demonstrated.  Shearing 
et  al.  [14]  conducted  volume  of  fluid  (VOF)  simulation  of  electro¬ 
chemical  reaction  and  diffusion  inside  Ni/YSZ  anode  obtained  from 
FIB-SEM  measurement.  In  the  above-mentioned  models,  however, 
each  of  the  solid  phases  is  assumed  to  be  either  perfect  electronic  or 
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Fig.  1.  (a)  Schematic  diagram  of  FIB-SEM  setting  and  (b)  reference  marks  milled  on  top  surface. 


ionic  conductors,  not  MIECs.  As  a  result,  electrochemical  reaction 
only  takes  place  at  the  TPB.  On  the  other  hand,  electrochemical 
reaction  can  be  active  at  wide  area  of  gas/solid  interface  in  the 
case  of  MIEC  electrodes,  since  both  electrons  and  ions  can  exist  in 
the  single  solid  phase.  In  order  to  design  and  optimize  MIEC  elec¬ 
trodes,  it  is  very  effective  to  use  numerical  simulation  tools  which 
can  describe  species  diffusions  and  charge  transfer  inside  complex 
electrode  microstructure. 

In  the  present  study,  the  three-dimensional  microstructure  of 
La0.6Sro.4Coo.2Feo.803_5  (LSCF6428)  cathode  is  quantified  by  a  dual 
beam  FIB-SEM.  Marching  cube  method  is  used  to  quantify  sur¬ 
face  area  from  voxel  structures.  Assuming  surface  electron  transfer 
to  be  the  rate  limiting  step,  reaction  current  is  expressed  by  a 
Butler-Volmer  like  form  [15].  Transport  equations  for  electron, 
oxide  ion  and  gas  species  coupled  with  charge  transfer  at  gas/solid 
interface  in  the  three-dimensional  field  is  solved  by  the  lattice 
Boltzmann  method.  Predicted  overpotential  is  compared  with 
the  experimental  data  for  validation.  Finally,  three-dimensional 
chemical  potential  and  current  vector  distributions  inside  LSCF 
microstructure  are  presented. 

2.  Sample  preparation  and  3D  microstructure 
reconstruction 

In  the  present  study,  electrolyte  supported  button  cell  with 
LSCF  cathode  and  Ni/YSZ  anode  is  used  (Japan  Fine  Ceramics  Co., 


Ltd.).  Press  molded  8YSZ  electrolyte  disk  was  sintered  at  1500°C, 
followed  by  the  GDC10  interlayer  sintering  at  1500°C.  NiO-8YSZ 
(60:40  wt%)  anode  and  LSCF6428  cathode  were  sintered  at  1300 
and  1000°C,  respectively. 

Observation  and  quantification  of  the  three-dimensional 
microstructure  of  the  LSCF6428  cathode  are  facilitated  by 
FIB-SEM  (Carl  Zeiss,  NVision  40)  [7].  The  sample  was  infil¬ 
trated  with  epoxy  resin  (Marumoto  Struers  KI<)  under  vacuum 
so  that  the  pores  could  be  easily  distinguished  during  SEM 
observation.  Fig.  1  schematically  shows  a  typical  setting  for 
the  FIB-SEM  measurement.  In  this  study,  an  in-lens  sec¬ 
ondary  electron  detector  was  used  with  acceleration  voltage  of 
1.5  kV.  By  automatically  repeating  FIB  milling  and  SEM  imag¬ 
ing  sequence,  a  series  of  cross  sectional  images  necessary  for 
three-dimensional  structure  reconstruction  is  acquired.  Reference 
marks  are  created  on  the  carbon  deposition  layer  on  the  sam¬ 
ple  surface  to  assist  the  alignment  of  SEM  images,  as  shown  in 
Fig.  1(b). 

Separation  of  solid  and  pore  phases  is  carried  out  for  each  SEM 
image.  Epoxy  infiltration  enables  easy  binarization  of  the  phases  as 
shown  in  Fig.  2(a).  As  the  LBM  simulation  requires  cubic  meshing, 
the  resolution  of  each  cross  sectional  image  was  coarsened  from 
13.96  nm/pixel  to  59.18nm/pixel  which  is  the  FIB  milling  pitch  as 
shown  in  Fig.  2(b).  The  phase  which  has  larger  area  portion  was  cho¬ 
sen  as  the  representative  phase  for  the  coarser  mesh.  Fig.  3  shows 
the  reconstructed  three-dimensional  microstructure.  Original  SEM 


Binarization 


Fig.  2.  (a)  Phase  separation  and  (b)  mesh  alignment  for  voxel  reconstruction  (white:  LSCF,  black:  pore). 
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Fig.  3.  Reconstructed  3D  cathode  microstructure:  (a)  pore  phase  and  (b)  LSCF  phase. 


images  correspond  to  x-y  cross  sections,  and  FIB  milling  is  carried 
out  in  z-direction. 

If  gas/solid  interface  is  simply  defined  by  the  voxel  surface,  sur¬ 
face  area  will  be  considerably  overestimated.  In  the  present  study, 
marching  cube  method  [16]  is  used  to  calculate  the  surface  area. 
Marching  cube  method  considers  21  patterns  to  calculate  inter¬ 
face  area  from  the  phase  information  of  neighboring  8  voxels  as 
shown  in  Fig.  4.  Flowever,  it  is  reported  that  marching  cube  method 
still  overestimates  surface  area  by  approximately  8%  in  the  case  of 
sphere  [17].  Table  1  shows  the  surface  area  of  the  reconstructed 
microstructure  calculated  by  the  marching  cube  method  and  stere- 
ology  [18].  Surface  area  is  calculated  for  divided  regions  shown  in 
Fig.  5  to  check  variation  between  the  samples.  Surface  areas  from 
x-divided  andy-divided  samples  show  very  small  variations,  while 
quite  large  variation  is  found  between  z-divided  samples.  It  is  con¬ 
sidered  that  the  error  caused  from  the  image  alignment  procedure 
in  the  FIB  milling  direction  (z-direction)  is  the  possible  reason  for 
the  variation  between  z-divided  samples.  Marching  cube  method 
gives  fairly  good  agreement  with  stereology  at  the  same  resolution 
of  59.18  nm/pixel.  Flowever,  stereology  results  show  dependency 
on  resolution,  and  it  should  be  considered  that  the  present  march¬ 
ing  cube  results  also  contain  grid  resolution  dependence. 


Table  1 

Surface  areas  calculated  by  marching  cube  method  and  stereology  for  different 
regions. 


Structure 

Volume  [fxm3] 

Surface  area 

density 

[|jim2/|jim3] 

Original 

23.731  x  10.593  x  13.493 

6.765 

x-axis  divided 

5.918  x  10.593  x  13.493 

6.641 

6.813 

6.866 

6.667 

y-axis  divided 

23.731  x  5.267  x  13.493 

6.666 

6.831 

z-axis  divided 

23.731  x  10.593  x  6.747 

7.145 

6.358 

Stereology  [18] 

(Resolution:  13.96  nm/px) 

7.734  (a:  0.859) 

(Resolution:  59.18  nm/px) 

6.035  (or;  0.536) 

(a) 


No  modified 
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Fig.  4.  Marching  cube  method  for  surface  area  calculation.  Black  and  white  spheres  represent  phase  information  of  each  voxel  and  estimated  interface  is  colored  in  blue,  (a) 
Example  with  one  black  phase  and  seven  white  phases  in  neighboring  eight  voxels,  (b)  21  interface  patterns  used  in  marching  cube  calculation.  (For  interpretation  of  the 
references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  the  article.) 
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Fig.  5.  Divided  structures:  (a)  x-axis  divided,  (b)y-axis  divided,  and  (c)  z-axis  divided. 


Table  2 

Gas  diffusion  parameters  used  in  Eqs.  (4)-(7). 


Substance 

M  [kg/mol] 

f  [A] 

e/k  [I<] 

o2 

31.9988  x  10“3 

3.54 

88.0 

n2 

28.0314  x  10“3 

3.68 

91.5 

3.  Numerical  analysis  by  the  lattice  Boltzmann  method 

3.1.  Diffusion  equations 


In  the  pore  phase,  diffusions  of  oxygen  and  nitrogen  are  solved 
based  on  the  dusty  gas  model  (DGM)  [19].  Convection  is  neglected 
and  constant  total  pressure  is  assumed  in  the  present  study  for  the 
first  order  approximation. 


Ni  x^yjNi  -y»Ni  = 

A,k  ^  DtJ 
J  #  i 


(1) 


where  y*  is  the  molar  fraction  and  Nz  is  the  molar  flux.  Subscripts  i 
and  j  present  gas  species,  such  as  oxygen  and  nitrogen.  For  constant 
total  pressure,  Graham’s  law  holds: 


Y)iiVMi=0-  (2) 


Therefore,  the  diffusion  equation  of  oxygen  can  be  expressed  as: 

-l  \ 


1  -  ayo2 

do2,n2 


+ 


1 


D02,k 


VG 


o2 


1  . 

:  ~4f Ireac’ 


(3) 


a  =  1  - 


Mq2 

Mn2 


(4) 


In  Eq.  (3),  Do2,n2  and  Do2,i<  rePresent  binary  and  Knudsen  diffu¬ 
sion  coefficients,  respectively: 


jo2,n2 


=  0.018833 


1  1 

+  ■ 


7'3/2 


Mq2  '  MN2  P^dG2 


o2,n2 


(5) 


„  2  /  8 RT 

°2’k  - 1  y 


where  is  the  collision  integral  given  as: 

/  Tk\  _01814 
^D  =  1.1336fyJ 


(6) 


(7) 


When  calculating  binary  diffusion  coefficient,  an  intermolecular 
force  constants  £  is  taken  as  an  arithmetic  mean  of  £q2  and  £n2. 
Geometric  mean  of  Sq2  and  £N2  is  used  for  e.  The  gas  parameters 
are  shown  in  Table  2.  The  mean  pore  radius  is  calculated  by  the 
maximum  sphere  inscription  method  (MSI)  [20].  In  MSI,  a  sphere 


with  radius  r  =  h/2  is  initially  inserted  in  each  pore  voxel  of  size  h. 
Then,  the  sphere  radius  is  increased  until  the  sphere  contains  at 
least  one  solid  phase  voxel  inside.  This  radius  is  assigned  as  pore 
radius  for  every  voxels  that  are  included  in  the  sphere.  Local  pore 
radius  is  given  by  taking  the  maximum  of  the  assigned  radii  at  each 
voxel.  In  the  present  LSCF  sample,  the  mean  pore  radius  of  r  =  89  nm 
is  obtained  by  taking  the  average  of  local  radii. 

Electronic  and  ionic  transport  equations  are  solved  in  the  solid 
MIEC  phase: 

V  (^pVAe-)  =  -ireac,  (8) 

V(^VA02-)  =  W,  (9) 

where  jle- ,  /xo2- ,  cre-  and  oq2-  are  electrochemical  potentials  and 
conductivities  of  electron  and  oxide  ion,  respectively.  In  a  mixed 
conductor,  <re-  and  crQ2-  depend  on  local  oxygen  potential  and 
temperature.  In  the  present  study,  cre-  and  chemical  diffusion  coef¬ 
ficient  of  oxide  ion  D  are  fitted  from  the  experimental  data  of 
Bouwmeester  et  al.  [21  ]. 

logio  Oe-  =  -0.0237(log10po2)2  + 0.0034  log10po2 

+  4.8126  (1073.15K),  (10) 


logio  ^e-  =  -0.0222(log10po2)2- 0.0169  log10po2 

+  4.8056  (1023.15K),  (11) 


log10<7e-  =  -0.0095(log10po2)2  -  0.0011  •  logio Po2 

+  4.8152  (973. 15K),  (12) 


logio  Oe-  =  -0.008(log,oPo2)2  - 0.0024  log10po2 

+  4.8447  (923. 15K),  (13) 


logio  D  =  -O.1765(log10po2)2  - 0.2724  logioPo2 

-9.2256  (1073. 15K),  (14) 


logio  D  =  -0.1884(log1Qpo2)2  -  0.3243  •  logi0Po2 

-9.4969  (1023.15  K),  (15) 


logio  D  =  — O.1882(log10po2)2  —  0.2491  logi0Po2 
-9.7676  (973.15  K), 


(16) 
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Fig.  6.  Chemical  diffusion  coefficient  of  LSCF6428  from  Bouwmeester  et  al.  [21]. 
Solid  lines  represent  fitting  curves,  Eqs.  (14)-(17). 


log10D  =  -O.1252(log10po2)2  -  0.2051  •  log10po2 

-9.9554  (923.15  K).  (17) 

Assuming  prevailing  electronic  conductivity  of  the  perovskites, 
chemical  diffusion  coefficient  D  is  related  to  the  ionic  conductivity 
<?0  2—  as: 


Local  oxygen  partial  pressure  in  the  solid  phase  is  calculated 
assuming  local  equilibrium: 

Mo  =  £02-  -2£e-,  (20) 

Mo  =  ^RTlogpo2.  (21) 

3.2.  Electrochemical  reaction  at  gas/solid  interface 

In  MIEC  electrodes,  charge  transfer  at  the  gas/solid  interface 
becomes  dominant.  According  to  Fleig  et  al.  [1 5,22],  electron  trans¬ 
fer  to  an  adsorbed  Oad  atom  at  the  gas/solid  interface  is  assumed  to 
be  the  rate  limiting  step  in  the  LSCF  cathode  reaction: 

Oad  +  e  ^  Oad.  (22) 

If  overpotential  is  applied  to  the  electrode,  gas/solid  interface 
potential  step  x  deviates  from  its  equilibrium  valuexeq  as: 

Ax  =  X  -  Xeq  =  <Pode  -  <Pad  -  (^ode  “  <0'  (23) 

Assuming  that  subsequent  elementary  steps  are  at  equilibrium 
and  ion  conduction  in  the  electrolyte  is  fast,  surface  potential  step 
A/  and  electrode  activation  overpotential  pact  can  be  related  as: 

RT  (  1  “  C 

2pact  =  Ax+T.n  T  _«  ^ 

\  ad  °ad 


1  a02-ae-  dn0  _  RTVm  91npo2 

4 F2  oq2-  +  ae-  dc0  *  8 F2  °°2~  dS 


(18) 


where  c0  is  the  concentration  of  oxygen  in  the  perovskite,  and 
Vm  =  35.5  x  10-6  m3/mol  is  the  molar  volume  of  LSCF.  Oxygen  non¬ 
stoichiometry  8  and  oxygen  partial  pressure  po2  is  again  fitted  from 
the  data  of  Bouwmeester  et  al.  [21]: 


— - =  -3.36260  x  10"5  ■  T  +  2.59403  x  10“2. 

9  In  po2 


(19) 


Fig.  6  shows  the  chemical  diffusion  coefficients  of  Bouwmeester 
et  al.  [21].  Solid  lines  in  Fig.  6  represent  the  fitting  curves  using 
Eqs.  (14)-(17).  As  described  in  Eq.  (19),  almost  linear  relationship 
is  found  between  8  and  logpo2  *  which  means  that  D  and  aQ2-  show 
nearly  the  same  trend  against  po2.  Profound  decline  in  chemical 
diffusion  coefficient  below  p02  of  10-2  bar  shown  in  Fig.  6  is  con¬ 
sidered  to  be  due  to  the  ordering  of  oxygen  vacancies,  which  is 
induced  by  the  increased  vacancy  concentration  and  the  enhanced 
interactions  between  the  vacancies  [21  ]. 


where,  O0-  denotes  surface  coverage  of  Oad.  If  the  surface  has 

high  number  of  surface  sites  and  surface  coverage  is  moderate, 
i.e.  O0-  «  i  ,  the  second  term  in  the  R.H.S  of  Eq.  (24)  can  be 

ad  ad 

neglected.  Then,  2pact  =  Ax  holds,  and  Butler-Volmer  like  equation 
with  an  additional  factor  of  two  in  the  exponents  can  be  used  for 
the  reaction  current  at  the  gas/solid  interface  [15]. 

'reac,2PB  =  io^2PB  jexp  (  ff'Uct'j  ~  exp  (-  ^^“^act)  )  •  (25) 

According  to  the  experimental  data  of  dense  LSCF  cathode  [22], 
Tafel  lines  with  slopes  of  1 .2 F/RT  and  1 .0 F/RT  in  anodic  and  cathodic 
regimes  are  observed.  In  the  present  study,  following  equation  is 
used  for  the  reaction  current  at  the  gas/solid  interface  as  a  function 
of  activation  overpotential  pac t: 

)  }  •  (26) 

Oxygen  partial  pressure  dependence  and  activation  energy  of 
the  lineal  exchange  current  z'o  is  fitted  from  the  experimental  data 


f  (\.2F  \  (  1  .OF 

*reac,2PB  =  Z0^2PB  |  exp  (  -j^"*7act  )  “  exp  f  --^pact 
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Fig.  7.  Total  cathode  overpotential. 
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Electrolyte 

(thickness  1.243  Tumi) 


Current  collector 
(thickness 


Fig.  8.  Computational  domain  (red:  current  collector;  gray:  electrolyte;  green:  LSCF).  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred 
to  the  web  version  of  the  article.) 


of  Esquirol  et  al.  [23]. 

i0  =  1.47  x  10 6  p°2* 2  •  exp  {~^f^)  •  (27) 

The  coefficient  1.47  x  106  in  Eq.  (27)  is  chosen  so  that  the  pre¬ 
diction  fits  most  to  the  experimental  data  shown  in  Section  4. 

Local  activation  overpotential  r]ac t  at  gas/solid  interface  is 
defined  as: 

Pact  =  ~2p  (2Ae-,MIEC  “  A02 *-,MIEC  +  2^^°^°2’Sas)  '  (28) 

The  total  cathode  overpotential  ^/cathode  is  obtained  by  subtract¬ 
ing  ohmic  losses  of  current  collector  (CC),  electrolyte  and  reference 
electrode  (RE)  from  the  total  overpotential,  defined  as  the  differ¬ 
ence  between  EMF  and  terminal  voltage  as  shown  in  Fig.  7. 


Table  3 

Numerical  conditions  for  the  present  work. 


Properties 

Value 

Operating  temperature,  T  [I<] 

973.15 

1023.15 

1073.15 

Total  pressure,  P  [Pa] 

1.013  x  105 

Gas  composition  (02:N2)  [mol%] 

100:0 

50:50 

20:80 

Current  density,  i  [A/cm2  ] 

0.01 

0.05 

0.1 

0.2 

^cathode  ~  ^Nernst  (^CC/S  ^RE/s)  ^7ohm,CC  ^7ohm,lyte  *7ohm,RE 

RT  f  Po2,cc\  ~  >. 

=  4Fl0g  [  p^ )  ~  F^'-’ RE/S  -  ^-,CC/S) 

~p  (Ae-,CC/S  -  Ae~ , cathode/ CC )  “  2p  (A02-,Cathode/lyte  _  A02-,lyte/RE) 
-^(Ae-,lyte/RE  “  Ae-,RE/s) 

RT  ( Po2,cc\  1  . 

~~  4 f l0g  l  p  R£  )  “  2f  l M'o2- , cathode /lyte  -  ^e-.cathode/CC  “  Mo,lyte/RE) 

=  ~2F  (^02~,cathode/lyte  _  2Ae-,cathode/CC  “  2  ^°SPo2,CC^  • 


(29) 


3.3.  Numerical  scheme 

The  lattice  Boltzmann  method  (LBM)  is  used  to  solve  Eqs.  (3),  (8) 
and  (9).  For  three-dimensional  LBM  simulation,  D3Q15  (1  =  1—15) 
or  D3Q19  (i  =  1—19)  models  are  commonly  used.  However,  in  the 
case  of  simple  diffusion  simulation,  it  has  been  shown  that  D3Q6 
(i  =  1  -6)  model  can  be  used  with  a  slight  loss  of  accuracy  [24].  Thus, 
the  D3Q6  model  is  used  in  this  study. 

The  lattice  Boltzmann  equation  with  the  LBGK  model  in  the 
collision  term  is  written  as  follows: 

/i(x  +  C,  At,  t  +  At)  =  fi(x,  t )  -  — L^[/i(x,  t)  _^(X,  t)]  +  w,  A t. 

(30) 


In  Eq.  {30),  fi  represents  the  density  distribution  function  of  gas, 
electron  or  ion  with  a  velocity  q  in  the  i-th  direction,  and  is  the 
Maxwellian  local  equilibrium  distribution, 

(3D 

1=1 


The  relaxation  time  t*  is  a  function  of  diffusion  coefficient  D, 
voxel  size  Ax  and  time  step  At  and  it  given  as: 


t*(x,  t)  =  0.5  + 


D(x,  t)At 

3  Ax2 


(32) 


In  this  study,  diffusion  coefficient  of  gas,  electronic  conductivity 
and  ionic  conductivity  are  spatially  non-uniform,  so  the  relaxation 
time  is  varied  locally  during  the  calculation. 

For  reducing  computational  time,  the  volume  is  halved  iny  and  z 
directions.  The  electrolyte  layer  is  added  atx  =  0  p,m,  and  the  current 
collector  (CC)  is  set  at  x  =  23.731  p-m.  Thicknesses  of  electrolyte  and 
CC  layers  are  1.243  pm  and  0.592  pm,  respectively.  The  resultant 
computational  domain  is  represented  by  432  x  90  x  1 15  voxels,  as 
shown  in  Fig.  8.  The  differences  of  specific  surface  area  and  vol¬ 
ume  fraction  between  the  reduced  size  structure  and  the  original 
structure  were  less  than  5%. 

Adiabatic  boundary  condition  is  applied  at  the  boundaries  of 
y  =  0,  5.267  pm  and  z  =  0,  6.747  pm.  Constant  gas  composition 
(Dirichlet  boundary)  is  given  at  the  current  collector  surface. 
Constant  electronic  and  ionic  current  flux  conditions  (Neumann 
boundary)  are  imposed  on  the  current  collector  and  electrolyte 
boundaries,  respectively.  A  zero-flux  boundary  is  imposed  on 
the  solid  surface  in  the  porous  media  by  applying  the  halfway 
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Fig.  9.  Comparison  of  cathode  overpotential,  (a)  Temperature  dependence  at  O2  =  20%,  and  (b)  oxygen  potential  dependence  at  T=  1023  K. 


bounceback  scheme  with  a  second-order  accuracy.  The  numerical 
conditions  are  shown  in  Table  3. 

4.  Computational  results 

Predicted  overpotential  is  compared  with  the  experimental 
results  in  Fig.  9.  Fig.  9(a)  shows  temperature  dependence  at 
02  =  20%,  and  (b)  shows  oxygen  potential  dependence  at  T=  1 023  K. 


The  error  bars  represent  standard  deviations  of  the  experi¬ 
mental  data.  Prediction  shows  quantitative  agreement  with  the 
experimental  data.  However,  the  predicted  overpotential  value 
at  O2  =  20%,  T  =  973  K,  i  =  0.2  A/cm2  case  was  considerably  large, 
approximately  0.19  V.  Since  this  value  is  far  beyond  the  vertical 
axis  range  of  Fig.  9(a),  the  predicted  point  is  not  shown  in  the  figure. 
This  overprediction  can  be  attributed  to  the  reduction  of  chemical 
diffusion  coefficient  of  LSCF  at  low  oxygen  partial  pressures  and 


Current  collector 


(b)  x  =  0.976  pm  (c)  x  =  5 .00 1  pm  (d)  jc  =  10.03 1  pm 


Fig.  10.  Oxygen  chemical  potential  /zo  distribution  in  LSCF  cathode  at  T=  1023  K,  O2  =  20%,  i  =  0.1  A/cm2. 
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Fig.  11.  Three-dimensional  current  stream  lines  in  LSCF  cathode  (red:  electronic  current;  blue:  ionic  current)  at  T=  1023  K,  02  =  20%,  i  =  0.1  A/cm2.  (For  interpretation  of  the 
references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  the  article.) 


low  temperatures  shown  in  Fig.  6.  It  is  considered  that  potential 
drop  near  the  electrolyte  side  of  the  LSCF  layer  lead  to  degradation 
of  ionic  conductivity,  which  further  decreased  the  potential.  This 
positive  feedback  caused  the  diverged  prediction  of  overpotential 
at  low  temperature  and  large  current  density  condition. 

Fig.  10  shows  three-dimensional  oxygen  chemical  potential  /x0 
distribution  at  T=1023K,  O2  =  20%,  i  =  0.1  A/cm2.  Note  that  color 


scales  in  Fig.  10(b)-(d)  for  the  cross  sectional  planes  are  magnified 
so  that  very  small  potential  difference  can  be  identified.  Oxygen 
chemical  potential  is  nearly  uniform  in  the  cross  sectional  plane, 
and  potential  gradually  drops  from  the  current  collector  to  the 
electrolyte.  Fig.  1 1  shows  three-dimensional  current  stream  lines. 
Electronic  and  ionic  current  stream  lines  are  nearly  straight  and 
parallel,  which  indicates  good  ionic  and  electronic  conductions  in 
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Fig.  12.  Chemical  potential  distribution  in  LSCF  cathode,  (a)  Dependences  on  current  density  at  02  =  20%  and  T=1023K,  (b)  dependence  on  temperature  at  02  =  20%  and 
i  =  0.1  A/cm2,  and  (c)  dependence  on  02  concentration  at  T=  1023  K  and  i  =  0.1  A/cm2. 
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Fig.  13.  Current  density  profiles  along  LSCF  thickness  direction.  Dependences  on  (a)  temperature,  and  (b)  bulk  O2  concentration. 


the  LSCF  phase.  Electronic  and  ionic  currents  reach  to  the  elec¬ 
trolyte  and  current  collector  sides,  respectively.  This  indicates  that 
the  electrochemical  reaction  is  taking  place  in  the  wide  electrode 
region. 

Fig.  12  shows  averaged  oxygen  chemical  potential  /x0  pro¬ 
file  along  the  LSCF  thickness  direction.  Oxygen  chemical  potential 
decreases  from  the  bulk  value  with  the  increase  of  current  den¬ 
sity  as  shown  in  Fig.  12(a).  Chemical  potential  decreases  mostly  in 
the  region  within  15  [xm  from  the  electrolyte.  The  bars  represent 
standard  deviation  (±1  sigma)  of  /x0  in  the  cross  sectional  plane. 
Thus,  very  uniform  potential  field  is  formed  in  the  cross  sectional 
plane.  Fig.  12(b)  shows  oxygen  chemical  potentials  at  i  =  0.1  A/cm2. 
As  shown  in  Fig.  12(b),  oxygen  chemical  potential  decreases  near 
the  electrolyte  as  temperature  decreases.  Since  bulk  /z0  is  higher  at 
lower  temperatures  due  to  the  temperature  dependence  of  /x0  as 
shown  in  Eq.  (21),  chemical  potential  profiles  incidentally  inter¬ 
sect  at  approximately  10|xm  from  the  electrolyte.  The  slope  of 
fio  is  steeper  for  100%  02  case  than  for  the  diluted  cases  at  the 
vicinity  of  the  electrolyte  as  shown  in  Fig.  12(c).  This  is  because 
exchange  current  increases  at  higher  po2  as  shown  in  Eq.  (27), 
while  ionic  conductivity  remain  nearly  unchanged  for  this  condi¬ 
tion  (po2  >  0.075  bar,  7=1023 1<). 

Due  to  the  oxygen  partial  pressure  dependence  of  diffusion 
coefficient  shown  in  Fig.  6,  decrease  in  oxygen  chemical  poten¬ 
tial  results  in  degradation  of  ionic  conductivity  near  the  electrolyte 
at  low  temperature  and  high  current  density  conditions.  This  is 
considered  to  be  the  main  reason  for  the  diverged  prediction  at 
02  =  20%,  7=973  K  and  i  =  0.2  A/cm2  case  as  discussed  above.  Fur¬ 
ther  investigations  on  the  p02  dependence  of  ionic  conductivity, 
especially  at  lowp02  and  low  temperatures,  should  be  required  for 
improving  the  accuracy  of  the  numerical  simulation  especially  for 
intermediate  temperature  SOFCs. 

Fig.  13  shows  ionic  and  electronic  current  distributions  along 
LSCF  thickness  direction.  As  can  be  seen  in  Fig.  13(a),  reactive 


region  becomes  thinner  as  temperature  decreases.  This  is  because 
of  the  degraded  ionic  conductivity  at  low  temperature.  On  the  other 
hand,  reactive  region  gets  thinner  as  bulk  oxygen  concentration  is 
increased,  as  shown  in  Fig.  13(b).  This  is  mainly  due  to  the  p02 
dependence  of  exchange  current  as  described  in  Eq.  (27).  Rela¬ 
tive  enhancement  of  surface  reaction  to  ionic  conduction  at  larger 
oxygen  concentration  reduces  reactive  region  thickness. 

5.  Conclusions 

Three-dimensional  microstructure  of  La0.6Sro.4Coo.2Feo.803_5 
(LSCF6428)  cathode  is  reconstructed  by  FIB-SEM,  and  its  over¬ 
potential  is  predicted  by  the  lattice  Boltzmann  method.  Surface 
area  estimated  by  the  marching  cube  method  showed  fairly  good 
agreement  with  those  from  stereology.  Electron  transfer  to  the 
adsorbed  oxygen  atom  on  the  gas/solid  interface  is  assumed  to  be 
the  rate  limiting  step  in  the  electrochemical  modeling.  The  cath¬ 
ode  overpotential  prediction  agreed  well  with  the  experimental 
data.  However,  predicted  overpotential  was  very  large  at  02  =  20%, 
7 =  973 1<  and  i  =  0.2  A/cm2  case.  This  can  be  attributed  to  the  decline 
in  oxygen  ionic  conductivity  at  low  po2.  Further  investigation 
on  oxygen  partial  pressure  dependence  of  ionic  conductivity  is 
required  for  improving  the  accuracy  of  the  numerical  simulation. 
Three-dimensional  chemical  potential  and  current  distributions 
are  uniform  and  smooth,  which  indicates  good  ionic  and  electronic 
conductions  as  well  as  wide  electrochemically  active  areas  inside 
the  LSCF  microstructure.  Present  method  will  be  an  effective  tool 
for  investigating  local  oxygen  potential  field  which  affects  local 
reactions,  diffusions  and  physical  properties  of  the  MIEC  cathodes. 
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